Physical origin of satellites in photoemission of doped graphene: an ab initio GW 

plus cumulant study 
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We calculate the photoemission spectra of suspended and epitaxial doped graphene using an 
ab initio cumulant expansion of the Green's function based on the GW self-energy. Our results 
are compared to experiment and to standard GW calculations. For doped graphene on a silicon 
carbide substrate, we find, in contrast to earlier calculations, that the spectral function from GW 
only does not reproduce experimental satellite properties. However, ab initio GW plus cumulant 
theory combined with an accurate description of the substrate screening results in good agreement 
with experiment, but gives no plasmaron (i.e., no extra well-defined excitation satisfying Dyson's 
equation). 



Introduction. — The isolation of graphene fH , an atomi- 
cally thin layer of carbon, in 2004 has sparked much inter- 
est and an enormous amount of scientific work. Among 
the most fascinating properties of this extraordinary ma- 
terial is its electronic structure: the low energy valence 
and conduction bands form two Dirac cones in the Bril- 
louin zone causing electrons in graphene to behave like 
massless relativistic particles with chiral character. 

Recently, several experimental groups probed the low- 
energy electronic structure of doped graphene on a sil- 
icon carbide (SiC) substrate using angle-resolved pho- 
toemission spectroscopy (ARPES) [2-4]. Surprisingly, 
the experimental spectra did not exhibit simple Dirac 
cones, but indicated a more complicated electronic struc- 
ture. There has been much discussion regarding the in- 
terpretation of the observed findings: while Lanzara and 
coworkers Q concluded that the SiC substrate induces 
a band gap in graphene, Rotenberg and coworkers 0, |1] 
argued that the observed electronic structure is intrinsic 
to graphene and is due to a plasmaron excitation caused 
by strong electron-carrier plasmon coupling. 

Both interpretations were backed by theoretical find- 
ings: Kim et al. Q carried out density- functional calcu- 
lations of graphene on a reconstructed SiC surface and 
found that the dangling bonds of the substrate bind 
to the graphene layer breaking the symmetry of the 
graphene sheet and causing a band gap to open. On the 
other hand, Polini et al. [6| and Hwang and Das Sarma 
f?\ carried out Green's function calculations on the Dirac 
Hamiltonian (model system with two linear bands) using 
the GW approximation to the self-energy and a sim- 
plified description of the substrate. These calculations 
found strong low-energy plasmaronic bands which obey 
the experimentally observed scaling law as a function of 
doping, and Bostwick et al. reported agreement of theory 
with experiment after introducing an adjustable param- 
eter attributed to the SiC substrate screening [3]. 

The agreement of GW calculations and experiment for 
satellite properties in doped graphene is surprising since 



it is known that for bulk solids, such as silicon or potas- 
sium [9- 11], electron correlations beyond the GW approx- 
imation are needed to accurately describe plasmon satel- 
lites. For bulk solids, the GW plus cumulant (GW-I-C) 
approximation 13, U], which includes significant vertex 
corrections, provides an accurate description of plasmon 
satellites, but does not find plasmarons. However, no ap- 
plication of the GW-I-C theory to a two dimensional sys- 
tem, such as doped graphene with its characteristic low- 
energy carrier plasmon, has been reported so far. Also, 
in contrast to bulk solids, the carrier plasmon dispersion 
in graphene is tunable and is both substrate and doping 
dependent. 

In this Letter, we investigate satellites in photoemis- 
sion spectra of doped suspended graphene and doped 
graphene on a SiC substrate using the ab initio cumu- 
lant expansion of the Green's function based on the GW 
approximation to the self-energy. We compare our find- 
ings with those from ab initio GW calculations and also 
with GW calculations based on the model Dirac Hamil- 
tonian. As a benchmark, we also calculated ab initio 
GW+C photoemission spectra for silicon. 

Methods. — The photoelectron current in an ARPES 
experiment with monochromatic photons of frequency i' 
and polarization is given by [12l ] 
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where k and a; are the momentum and binding energy 
of an electron in band n, /(w) is the Fermi-Dirac dis- 
tribution and /q includes the absorption cross section of 
the incident photons. Also, Ank[u) = l/7r)ImGjifc(a;)l 
denotes the spectral function with Gnk{^) being the in- 
teracting one-particle Green's function. 

Usually, Gnkiy) is obtained by solving Dyson's equa- 
tion G-l{u^) = G-X^{i^) - + with GoM^) 

and V^^ denoting a mean-field Green's function and 
exchange-correlation potential, respectively, and I]„fc(cj) 
is the self-energy. In this work, we employ the ab initio 
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GW approach to the self-energy [8|, [iSj . 

While describing quasiparticle properties in many ma- 
terials with high accuracy [l^l , the G W ap proximation is 



less reliable for satellite properties : for the spec- 



tral function of a core electron interacting with plasmons, 
GW predicts a single satellite instead of a satellite series 
with decreasing spectral weight and also greatly overesti- 
mates the binding energy of the satellite structures. The 
cumulant expansion ^Q, 14, [l^ of G„fe(a;) cures these 
deficiencies by including significant vertex corrections be- 
yond GW: it provides the exact solution for a core elec- 
tron interacting with plasmons In the cumulant ap- 
proach, the Green's function for a hole is expressed as 



(2) 



where e„fc denotes the mean-field orbital energy and 
Cnk{i) denotes the cumulant. This expression for the 
Green's function is obtained after the first iteration of 
the self-consistent solution of its equation of motion as- 
suming a simple quasiparticle form for the starting guess 

M- 

The cumulant can be separated into a quasiparticle 
part C'^f^{t) and a satellite part C!^{t) given formally in 
terms of the self-energy by (for t <Q) 



ImS„fc(w) 
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where fj, denotes the chemical potential, r/ is a positive 
infinitesimal, E^k = ^nk + ^nk^Enk) - V^{: is the quasi- 
particle energy and Ej|;^,(a;) is defined through the rela- 
tion 
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For a given level of approximation for E, the cumu- 
lant theory yields an improved Green's function through 
Eqs. dSKS]). In the present study, S is obtained from ah 
initio GW theory [13j which is known to describe quasi- 
particle properties in Si, G and related materials accu- 
rately thus providing a good starting point for the cumu- 
lant theory. 

Silicon. — To benchmark the accuracy of the ah initio 
GW-I-C approach, we first applied it to silicon where ac- 
curate angle-integrated photoemission data is available 
llj . We carried out one-shot full frequency GoWn cal- 
culations for E using the Berkeley GW package [16j [see 
Supplementary Materials (SM) for details] and then eval- 
uated the cumulant spectral function. 

The inset of Fig. [T] compares the spectral functions 
from ah initio GW and ah initio GW-I-C theory for the 
lowest valence band at the F point of silicon. While the 
GW spectral function exhibits a strong plasmaron peak 
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FIG. 1: Comparison of theoretical results and experiment 
datap^ for the angle-integrated photoemission spectrum of 
silicon. The inset shows the spectral function of the lowest 
valence band at the F-point (Fi„) in silicon from ab initio GW 
(dashed blue curve) and ab initio GW-I-C (solid red curve) 
theories. The energy is measured from the top of the valence 
band 



due to an additonal spurious solution of Dyson's equation 
[§-11 1 (see SM) separated by ~ 23 eV from the quasipar- 
ticle peak, the GW-f-C spectral function shows a more 
shallow peak separated by ^ 16 eV from the quasiparti- 
cle peak. This separation agrees well with the experimen- 
tal plasmon energy of 16.6 eV 17[ in silicon indicating a 
much weaker interaction between electrons and plasmons 
than predicted by GW. 

Figure [1] shows the angle-integrated photoemission 
spectrum from the ah initio GW-I-G theory and com- 
pares it to the ah initio GW result and experiment. We 
have included matrix element effects, the inelastic back- 
ground and extrinsic plasmon losses similar to the work 
of Guzzo and coworkers [HI (see SM). Our full frequency 
ah initio GW-I-G theory reproduces accurately the loca- 
tion, height and width of the first plasmon peak and it 
does not show a significant second peak. 

Doped Graphene. — In contrast to silicon where the 
plasmons are high-energy excitations, doped graphene 
exhibits a low-energy, dispersive plasmon branch due 
to added carriers from doping with a characteristic 2D 
square root of q dispersion relation 0, 0, 0] i the so-called 
carrier plasmon. Due to the special dispersion relations of 
the carrier plasmons and of the Dirac fermions, the self- 
energy of doped graphene exhibits very sharp features 
in (fc, w)-space necessitating extremely fine k-point 
and frequency grids for convergent results: our calcula- 
tion of T,nk{^) and Anki^) employs a full frequency sam- 
pling at 0.05 eV intervals and a 1440 x 1440 x 1 k-point 
grid (see SM). 

Figure [2l^a) shows the ah initio GW-I-G spectral func- 
tion at the Dirac point of isolated graphene with elec- 
tron doping resulting in a charge density n = —12 A x 



10 cm . Also shown are the ah initio GW spectral 
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function and the GW spectral function from the hnear- 
band model[^ (In the hnear-bands model, graphene 
is described by the Dirac Hamiltonian with a bandstruc- 
ture consisting of only two linear bands at the K and 
K' points of the Brillouin zone.) For the linear-band 
model, we included screening contributions from all other 
bands by adding the analytically parametrized cRPA di- 
electric function calculated within DFT by excluding the 
Pz bands, obtained by Wehling and coworkers (Eq. (3) 
of Ref.(l8l|). to the calculated dielectric function of the 
linear- band model [3, [i^l • We verified the accuracy of 
the cRPA dielectric function of Wehling et al. by com- 
paring it to the additional background dielectric constant 
K(q) necessary to match the carrier plasmon peak of the 
imaginary part of the interacting susceptibility from the 
linear-band model to the corresponding peak in the ab 
initio calculation [see inset of Figl2{b)]. 

Despite the simplification of the graphene band struc- 
ture in the linear-band model, the resulting GW spectral 
function is similar to the ab initio GW result provided 
that the above K{q) is used: both theories give a strong 
quasiparticle peak and a smaller satellite peak. The sepa- 
ration of the peaks is similar, 0.83 eV in ab initio GW and 
0.89 eV in the linear-band model. However, the shapes 
of the features differ noticeably: in the ab initio GW 
theory, the quasiparticle peak is much broader than in 
the linear-band model, while the satellite peak is sharper 
revealing the importance of a realistic bandstructure for 
high doping levels. The GW approximation leads to a 
spurious plasmaron solution of Dyson's equation giving 
rise to the peak near -2 eV in Fig. 2(a) (see SM). The in- 
clusion of higher order electron interaction effects in the 
ab initio GW-I-C calculations show important changes in 
the spectral function in Fig. [2j i) there is a significant 
shift of the amplitude from the quasiparticle peak to the 
satellite peak, and ii) the peak separation is reduced to 
0.62 eV, a change of ~ 30 percent. The inset in Fig- 
me^a,) shows that the ab initio GW-I-C spectral func- 
tion also exhibits an additional shoulder-like structure at 
cj ^ —2.3 eV resulting from the coupling of the hole to 
two plasmons. So far no experimental ARPES spectrum 
with sufficient accuracy to observe satellite features has 
yet been reported for suspended graphene. 

Bostwick and coworkers had carried out ARPES ex- 
periments on hydrogen-intercalated epitaxial graphene 
grown on SiC(OOOl) 0, [22|. To include the screen- 
ing effect of the substrate we model the substrate by 
a 10-atomic-layer thick slab of hydrogen-terminated 4H- 
SiC(OOOl) employing the V3 x V^i?30° SiC surface unit 
cell [23l425| [see inset of Figure [31[b)] resulting in a su- 
percell geometry with neighboring graphene sheets sep- 
arated by 54 a.u.. To determine the SiC-graphene dis- 
tance we carried out density-functional theory calcula- 
tions with the empirical van der Waals correction of 
Grimme [l^ and obtain a separation of 3.86 A which 
agrees well with the findings of Soltys and coworkers [2^ . 



(a) 



linear-band model GW 
ab initio GW 
ab initio GW+G 




(b) 



E 



0.008 
0.007 
0.006 
0.005 
0.004 
0.003 
0.002 
0.001 




ab initio isolated graphene 
ab initio graphene on SiC 
linear-band model with K(q■^)=2.2 




a (eV) 



FIG. 2: (a): Comparison of ab initio GW-I-C, ab initio GW, 
and linear-band model GW spectral functions at the Dirac 
point for isolated graphene with electron doping correspond- 
ing to a charge density n — —VIA x 10^'^ cm~^. Energies 
are measured relative to the Fermi energy. The inset shows 
an additional shoulder-like feature (indicated by arrow) in 
the ah initio GW-I-C spectral function, (b): Imaginary part 
of the head of the interacting susceptibility (in a supercell 
calculation Imxoo(q,<^) ~ ^^^oo il' ^) /'"(l) /''^a with v{q) de- 
noting the truncated Coulomb interaction and Ug the density 
of graphene sheets per unit length) for isolated graphene and 
for graphene on a SiC substrate. Note that ImYoo(q,aj) is 
proportional to the dynamical structure factor S{q,Lo) [211 ]. 
We also show the result from the linear-band model with 
a background dielectric of K{qi) — 2.2. The results are for 
n = —5.9 X 10^^ cm"^ and gi = 0.022 a.u.. The inset shows 
the background dielectric function n{q) used in the strictly 
2D linear-band model for isolated graphene (solid circles), 
graphene on SiC (solid squares) and the cRPA result for iso- 
lated graphene obtained by Wehling et al. (solid line') [l8l| . 



Because of the passivation of the surface dangling bonds 
by hydrogen atoms, we also do not find a surface induced 
band gap opening in the graphene sheet [25 1. 

We then carried out GW calculations on the 
graphene-f substrate system using a truncated Coulomb 
interaction fl6|. Due to the immense size of the system 
under consideration, we separately calculate the contri- 
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butions to the polarizability matrix from the graphene 
sheet and the substrate before adding them up (see SM). 

Figure Elb) shows the imaginary part of the interact- 
ing susceptibihty of graphene on H-passivated SiC for 
n = —5.9 X 10"'^^ cm~^ at a small wave vector qi = 
0.022 a.u. and compares it to the isolated graphene re- 
sult. The substrate screening leads to a reduction of the 
carrier plasmon frequency (the position of the sharp peak 
at ^ 0.7 eV). Figure [2l^b) also shows that the plasmon 
peak of the linear-band model coincides with the ab ini- 
tio result if the 2D bare Coulomb interaction is divided 
by a background dielectric constant K{qi) — 2.2 which 
now contains both the intrinsic graphene screening pro- 
cesses discussed earlier and screening contributions from 
the substrate. This deduced ab initio value of k is much 
smaller than the value used by Bostwick and coworkers 
Q who treated k as a q-independent fitting parameter 
and found that k ~ 5.15 (using the LDA Fermi velocity 
to convert the coupling strength ac = 0.5 to a back- 
ground dielectric constant) gives the best description of 
the experimental spectra within the GW approximation 
only. We find that k increases with q, but remains much 
smaller than the value needed by Bostwick and coworkers 
to fit experiment [see inset of Fig. [2fb)]. 

Figure |3lja) compares ab initio GW+C and GW spec- 
tral functions for graphene on SiC at the Dirac point 
with experiment for n = —5.9 x 10^'^ cm~^. All spec- 
tral functions exhibit a quasiparticle peak and a satel- 
lite: their separation is 0.44 eV in GW and 0.27 eV in 
the GW-I-C theory. The ab initio GW+C theory agrees 
very well with the experimental separation of 0.30 eV W\. 
Also, the shape of the GW-I-C spectral function resembles 
the experimental spectrum more closely than the GW re- 
sult. The inset of Fig. ^a) shows that GW predicts an 
additional plasmaron solution to Dyson's equation. In 
contrast, no plasmaron solution is found if the vertex- 
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corrected self-energy ST^gw+c = Gq — G, 



GW+C 



is used. 



We attribute the remaining difference between the exper- 
imental spectrum and the computed spectral function to 
extrinsic losses (i.e., emission of plasmons by the pho- 
toexcited electron), which are not included in our theory 
(see SM for a preliminary treatment of extrinsic losses). 

Bostwick and coworkers have shown that the exper- 
imental spectral functions exhibit _an important scaling 
behavior as function of doping [3]: when the exper- 
imental quasiparticle-satellite separation is divided by 
the Fermi level Ep the resulting number is independent 
of doping for a wide range of doping levels, a conse- 
quence of the linear dispersion of electrons in graphene. 
To check whether the ab initio GW-I-C method repro- 
duces this scaling behavior we carried out calculations 
for n = —1.5 X 10^'^ cm~^. The results are shown in 
Fig. El^b). Again, we find that the separation of quasi- 
particle and satellite peaks is described accurately by the 
ab initio GW-I-C method. 

In conclusion, we have carried out a first-principles 
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FIG. 3: Comparison of ab initio GW-I-C and GW spectral 
function (in 1 /eV) for graphene on SiC at the Dirac point with 
experiment (experimental spectra are given in counts/second 
in arbitrary units) for (a) n = —5.9 x lO'^^ cm~^ and (b) 
n = —1.5 X 10^^ cm"'^. The arrows indicate the satellite 
position in the theoretical curves. The inset of (a) shows 
Re5E(a;) - to + slda with = Gq^ - for GW (blue 
curve) and GW-I-C (red curve) as a function of oj. The arrow 
denotes the plasmaron solution in GW, which does not exist in 
GW-I-C. The inset of (b) shows the geometry of graphene on 
a hydrogen-terminated 4H-SiC (0001) surface; carbon atoms 
are brown, silicon atoms are blue and hydrogens are white. 



study of doped graphene on hydrogen- intercalated SiC, 
going beyong standard GW calculations, which explains 
the experimental findings of Bostwick and coworkers Q . 
The work reproduces the experimental satellite proper- 
ties without finding a plasmaron showing the importance 
of an advanced treatment of electron correlations via the 
ab initio GW-I-C method and an accurate modelling of 
substrate screening. We do not find a substrate induced 
band gap or a plasmaron solution. We also explain the 
previously reported good fit of GW calculations to ex- 

geriment using an simplified modelling of the substrate 
|. In these calculations, the lack of sufficient electron 
correlations at the GW level leads to an extra solution to 
Dyson's equation which gives rise to an overestimation of 
the quasiparticle-satellite separation. However, these cal- 
culations also overestimate the effect of substrate screen- 
ing resulting in spurious agreement with experiment. 
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